---
title: "`r config.population.name` stock assessment on data-limited CMSY model (`r year.start` - `r year.terminal`) in `r config.population.area`"
author: "Developer: Piatinskii M., Report build by: `r config.report.author`"
date: 'Report build date: `r Sys.time()`'
output: 
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## 1. Model information
**CMSY** - The CMSY model developed by Froese et al. 2017 employs a stock reduction analysis using priors for r based on resilience, K based on maximum catch and the r priors, and start, intermediate, and final year saturation based on a set of simple rules. It also allows users to revise the default priors based on expert knowledge. The SRA employs a Schaefer biomass dynamics model and an algorithm for identifying feasible parameter combinations to estimate biomass, fishing mortality, and stock status (i.e., B/BMSY, F/FMSY) time series and biological/management quantities (i.e., r, K, MSY, BMSY, FMSY).

To perform CMSY model only catch information required. 

Reference: [Froese R, Demirel N, Coro G, Kleisner KM, Winker H (2017) Estimating fisheries reference points from catch and resilience. Fish and Fisheries 18(3): 506-526.](http://onlinelibrary.wiley.com/doi/10.1111/faf.12190/abstract)


## 2. Input data
This section present the input data for stock assessment procedure. 

Required input data:

  - C - catch information over years. Skips not allowed.

```{r}
knitr::kable(data, caption = "Table 2.1. Input data", align = "l")
```

## 3. Model tuning
Model resilience tuning: **`r config.population.resilience`**

CMSY model performed using next cmd: 

```{r eval = FALSE, echo = TRUE}
cmsy <- cmsy2(year=data$year, catch=data$catch, resilience = config.population.resilience)
```

## 4. Results
There summary CMSY modelling results shown. Feel free to use it.

### 4.1 Estimates - B, F

```{r echo = FALSE}
res <- cmsy$ref_ts %>%
  select(year, b, b_lo, b_hi, f, f_lo, f_hi) %>%
  mutate_at(c(3,4), ~round(., 0)) %>%
  mutate_at(c(6,7), ~round(., 3)) %>%
  unite("b.ci95", b_lo:b_hi, sep = " - ", remove = TRUE) %>%
  unite("f.ci95", f_lo:f_hi, sep=" - ", remove = TRUE)

knitr::kable(res, caption = "Table 4.1.1. Biomass and fishing mortality estimates", align="l")
```

Column description:

  - b - biomass estimation
  - b.ci95 - biomass confidence interval at p = 0.95 level
  - f - fishing mortality estimation
  - f.ci95 - fishing mortality conf.interval at p = 0.95
  
```{r echo = FALSE, fig.cap="Fig. 4.1.1. B estimates with 95% confidence interval"}
ggplot(data = cmsy$ref_ts, aes(x = year)) + 
  ylim(c(0, max(cmsy$ref_ts$b_hi))) + 
  geom_line(aes(y = b)) + 
  geom_point(aes(y = b)) + 
  geom_line(aes(y = b_lo), linetype = "dashed", alpha = 0.15) + 
  geom_line(aes(y = b_hi), linetype = "dashed", alpha = 0.15) + 
  geom_ribbon(data = cmsy$ref_ts, aes(ymin = b_lo, ymax = b_hi), alpha = 0.15) + 
  xlab("Year") + 
  ylab("Biomass (B)") + 
  theme_minimal()
```

```{r echo = FALSE, fig.cap="Fig. 4.1.2. F estimates with 95% confidence interval"}
ggplot(data = cmsy$ref_ts, aes(x = year)) + 
  ylim(c(0, max(cmsy$ref_ts$f_hi))) + 
  geom_line(aes(y = f)) + 
  geom_point(aes(y = f)) + 
  geom_line(aes(y = f_lo), linetype = "dashed", alpha = 0.15) + 
  geom_line(aes(y = f_hi), linetype = "dashed", alpha = 0.15) + 
  geom_ribbon(data = cmsy$ref_ts, aes(ymin = f_lo, ymax = f_hi), alpha = 0.15) + 
  xlab("Year") + 
  ylab("Fishing mortality (F)") + 
  theme_minimal()
```
  
### 4.2 Reference points
MSY reference point stategy done. Mean reference points shown in table 4.2.1.

```{r echo=FALSE}
  cmsy$ref_pts[-2:-1,] %>%
  mutate_if(is.numeric, ~round(., 3)) %>%
  mutate_if(is.numeric, ~format(., scientific = FALSE)) %>%
  knitr::kable(., caption="Table 4.2.1. MSY reference point estimates", align="l", row.names=FALSE, digits=4)
```

Column description:

  - param - reference point name column
  - est - estimated value
  - lo - lower confidence level
  - hi - upper (higher) confidence level

Remember, Bmsy calculated once (one time) in whole time seria. Fmsy point calculated for every one year but not presented here.

### 4.3 Estimates VS reference points

```{r echo = FALSE, fig.cap="Fig. 4.3.1. B/Bmsy proportion in time vector"}
ggplot(data = cmsy$ref_ts, aes(x = year)) + 
  ylim(c(0, max(cmsy$ref_ts$bbmsy_hi))) + 
  geom_line(aes(y = bbmsy)) + 
  geom_point(aes(y = bbmsy)) + 
  geom_line(aes(y = bbmsy_lo), linetype = "dashed", alpha = 0.15) + 
  geom_line(aes(y = bbmsy_hi), linetype = "dashed", alpha = 0.15) + 
  geom_ribbon(data = cmsy$ref_ts, aes(ymin = bbmsy_lo, ymax = bbmsy_hi), alpha = 0.15) + 
  geom_hline(yintercept = 1, linetype = "dashed", size = 0.3) + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = 1.05, label = "Btr", color = "darkgray") +
  geom_hline(yintercept = config.blim_frommsy, linetype = "dashed", size = 0.2) + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = config.blim_frommsy-0.05, label = "Blim", color = "darkgray") +
  xlab("Year") + 
  ylab("B/Bmsy reference") + 
  theme_minimal()
```

```{r echo = FALSE, fig.cap="Fig. 4.3.2. F/Fmsy proportion in time vector"}
ggplot(data = cmsy$ref_ts, aes(x = year)) + 
  ylim(c(0, max(cmsy$ref_ts$ffmsy_hi))) + 
  geom_line(aes(y = ffmsy)) + 
  geom_point(aes(y = ffmsy)) + 
  geom_line(aes(y = ffmsy_lo), linetype = "dashed", alpha = 0.15) + 
  geom_line(aes(y = ffmsy_hi), linetype = "dashed", alpha = 0.15) + 
  geom_ribbon(data = cmsy$ref_ts, aes(ymin = ffmsy_lo, ymax = ffmsy_hi), alpha = 0.15) + 
  geom_hline(yintercept = 1, linetype = "dashed", size = 0.3) + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = 1.05, label = "Ftr", color = "darkgray") +
  xlab("Year") + 
  ylab("F/Fmsy reference") + 
  theme_minimal()
```

```{r echo = FALSE, fig.cap="Fig. 4.3.3. Catch, biomass, MSY and Bmsy in time vector"}
ggplot(data = cmsy$ref_ts, aes(x = year)) + 
  ylim(0, max(cmsy$ref_ts$b_hi)) + 
  xlim(min(cmsy$ref_ts$year), max(cmsy$ref_ts$year)) + 
  geom_line(aes(y = b, color = "B & Bmsy")) + 
  geom_point(aes(y = b, color = "B & Bmsy")) + 
  geom_ribbon(data = cmsy$ref_ts, aes(ymin = b_lo, ymax = b_hi), fill = "blue", alpha = 0.05) + 
  geom_hline(yintercept = cmsy$ref_pts$est[5], linetype = "dotted", size = 0.3, color = "blue") + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = cmsy$ref_pts$est[5]*1.03, label = "Btr", color = "blue") +
  geom_hline(yintercept = cmsy$ref_pts$est[5]*config.blim_frommsy, linetype = "dashed", size = 0.3, color = "blue") + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = cmsy$ref_pts$est[5]*config.blim_frommsy*1.03, label = "Blim", color = "blue") +
  geom_line(aes(y = catch, color = "C & CMSY")) + 
  geom_point(aes(y = catch, color = "C & CMSY")) + 
  geom_hline(yintercept = cmsy$ref_pts$est[3], linetype = "dashed", size = 0.3, color = "red") + 
  geom_text(x = max(cmsy$ref_ts$year)-2, y = cmsy$ref_pts$est[3]*1.03, label = "MSY", color = "red") +
  labs(x = "Year", y = "Biomass / Catch", color = "") + 
  scale_color_manual(values = c("Blue", "Red")) +
  theme_minimal() + 
  theme(legend.position = c(0.09, 0.92))
```

### 4.4 Model r/K parametrization
```{r echo=FALSE}
  cmsy$ref_pts[1:2,] %>%
  mutate_if(is.numeric, ~round(., 3)) %>%
  mutate_if(is.numeric, ~format(., scientific = FALSE)) %>%
  knitr::kable(., caption="Table 4.4.1. Schaefer model r/K optimum found", align="l", row.names=FALSE, digits=4)
```

Column description:

  - param - param name column
  - est - estimated value
  - lo - lower confidence level
  - hi - upper (higher) confidence level

## 5. Forecast scenarious
This section contains forecast scenarious with different TAC levels. 

### 5.1. Scenarious data
Input scenarious data shown 
```{r echo=FALSE}
scenarios %>%
  knitr::kable(., caption = "Table 5.1.1. Forecast scenarious TAC (catch values)", align = "l", row.names = F, digits = 2)
```

### 5.2. Estimates - B, F
Biomass forecast features shown
```{r echo=FALSE, fig.cap="Fig. 5.2.1. Biomass forecast at TAC scenarious"}
sc <- names(forecast.b[,-1])
res.retro <- cmsy$ref_ts[(nrow(cmsy$ref_ts)-1):(nrow(cmsy$ref_ts)),]
df <- data.frame(year = c(res.retro$year, forecast.b$year))
for (s in sc) {
  df[,s] <- c(res.retro$b, forecast.b[,s])
}
p <- ggplot(data = df, aes(x = year)) + 
  ylim(0, max(df[,-1], cmsy$ref_pts[5, "est"])) + 
  geom_hline(yintercept = cmsy$ref_pts[5, "est"], linetype = "dotted") + 
  geom_hline(yintercept = cmsy$ref_pts[5, "est"] * config.blim_frommsy, linetype = "dashed") + 
  geom_text(x = max(df$year)-1, y = cmsy$ref_pts[5, "est"]*1.03, label = "Btr") + 
  geom_text(x = max(df$year)-1, y = cmsy$ref_pts[5, "est"]*0.5*1.03, label = "Blim")

for (s in sc) {
  p = p + geom_line(aes_string(y = s, color = shQuote(s)))
  p = p + geom_point(aes_string(y = s, color = shQuote(s)))
}
p = p + theme_minimal() + 
  labs(x = "Year", y = "Biomass (B)", color = "Scenario") + 
  theme(legend.position = c(0.03, 0.85))
plot(p)
```

```{r echo=F, fig.cap="Fig. 5.2.2. Fishing mortality scenarious"}
sc <- names(forecast.f[,-1])
res.retro <- cmsy$ref_ts[(nrow(cmsy$ref_ts)-1):(nrow(cmsy$ref_ts)),]
df <- data.frame(year = c(res.retro$year, forecast.b$year))
for (s in sc) {
  df[,s] <- c(res.retro$f, forecast.f[,s])
}
p <- ggplot(data = df, aes(x = year)) + 
  ylim(0, max(df[,-1], cmsy$ref_pts[4, "est"])) + 
  geom_hline(yintercept = cmsy$ref_pts[4, "est"], linetype = "dashed") + 
  geom_text(x = max(df$year)-1, y = cmsy$ref_pts[4, "est"]*1.03, label = "Ftr")

for (s in sc) {
  p = p + geom_line(aes_string(y = s, color = shQuote(s)))
  p = p + geom_point(aes_string(y = s, color = shQuote(s)))
}
p = p + theme_minimal() + 
  labs(x = "Year", y = "Fishing mortality (F)", color = "Scenario") + 
  theme(legend.position = c(0.03, 0.85))
plot(p)
```

Absolute biomass estimates: 
```{r echo=F}
if (config.forecast.use) {
  forecast.b %>%
    knitr::kable(., caption = "Table 5.2.1. Biomass forecast by scenarious", align = "l", row.names = F, digits = 2)
}
```

Fishing mortality estimates: 
```{r echo=F}
if (config.forecast.use) {
  forecast.f %>%
    knitr::kable(., caption = "Table 5.2.1. Biomass forecast by scenarious", align = "l", row.names = F, digits = 3)
}
```

### 5.3. Relative ref.pts
Relative biomass & fishing mortality to target levels MSY shown.
```{r echo=F, fig.cap="Fig. 5.3.1. Relative biomass B/Bmsy forecast scenarious"}
sc <- names(forecast.bbmsy[,-1])
res.retro <- cmsy$ref_ts[(nrow(cmsy$ref_ts)-1):(nrow(cmsy$ref_ts)),]
df <- data.frame(year = c(res.retro$year, forecast.bbmsy$year))
for (s in sc) {
  df[,s] <- c(res.retro$bbmsy, forecast.bbmsy[,s])
}
p <- ggplot(data = df, aes(x = year)) + 
  ylim(0, max(df[,-1], 1)) + 
  geom_hline(yintercept = 1, linetype = "dashed") + 
  geom_text(x = max(df$year)-1, y = 1.03, label = "Btr")

for (s in sc) {
  p = p + geom_line(aes_string(y = s, color = shQuote(s)))
  p = p + geom_point(aes_string(y = s, color = shQuote(s)))
}
p = p + theme_minimal() + 
  labs(x = "Year", y = "B/Bmsy", color = "Scenario") + 
  theme(legend.position = c(0.03, 0.85))
plot(p)
```

```{r echo=F, fig.cap="Fig. 5.3.2. Relative fishing mortality F/Fmsy forecast scenarious"}
sc <- names(forecast.ffmsy[,-1])
res.retro <- cmsy$ref_ts[(nrow(cmsy$ref_ts)-1):(nrow(cmsy$ref_ts)),]
df <- data.frame(year = c(res.retro$year, forecast.ffmsy$year))
for (s in sc) {
  df[,s] <- c(res.retro$ffmsy, forecast.ffmsy[,s])
}
p <- ggplot(data = df, aes(x = year)) + 
  ylim(0, max(df[,-1], 1)) + 
  geom_hline(yintercept = 1, linetype = "dashed") + 
  geom_text(x = max(df$year)-1, y = 1.03, label = "Ftr")

for (s in sc) {
  p = p + geom_line(aes_string(y = s, color = shQuote(s)))
  p = p + geom_point(aes_string(y = s, color = shQuote(s)))
}
p = p + theme_minimal() + 
  labs(x = "Year", y = "F/Fmsy", color = "Scenario") + 
  theme(legend.position = c(0.03, 0.85))
plot(p)
```

## 6. Diagnostics
Some avaiable model diagnostics are summarized here.

### 6.1. Retrospective analysis
There is summary retrospective figure presented. Retrospective procedure are syntetic continious reducing time seria by 1 year.

```{r echo=FALSE, fig.cap="Fig. 6.1.1. Retrospective B analysis"}
p <- ggplot(data = cmsy$ref_ts, aes(x = year)) +
  ylim(0, b.max) +
  geom_line(aes(y = b, color = "Bdefault"), size = 1.2)

for (i in 1:config.retro.years) {
  d <- as.data.frame(retro.data[[i]])
  col <- paste0("B", max(d$year))
  p = p + geom_line(data = d, aes_string(y = d$b, color = shQuote(col)))
}

p = p + 
  theme_minimal() + 
  labs(x = "Year", y = "Biomass (B)", color = "Retrospective") + 
  theme(legend.position = c(0.04, 0.85))
plot(p)

```


```{r echo=FALSE, fig.cap="Fig. 6.1.2. Retrospective F analysis"}
p <- ggplot(data = cmsy$ref_ts, aes(x = year)) +
  ylim(0, f.max) +
  geom_line(aes(y = f, color = "Fdefault"), size = 1.2)

for (i in 1:config.retro.years) {
  d <- as.data.frame(retro.data[[i]])
  col <- paste0("F", max(d$year))
  p = p + geom_line(data = d, aes_string(y = d$f, color = shQuote(col)))
}

p = p + 
  theme_minimal() + 
  labs(x = "Year", y = "Fishing mortality (F)", color = "Retrospective") + 
  theme(legend.position = c(0.04, 0.85))
plot(p)
```

### 6.2. Mohn-rho test
The basic ICES (2018) procedure to determine model stability over years - retrospective Mohn rho test. Mohn's rho test calculate relative bias for latest 5 years (default ICES procedure - 5 years) retrospective variations on scale -1 ... +1. Low negative values of *rho* leads to underestimate factor, high positive values - to overestimation (remark: values lowest -0.4 or higher +0.4 shows high variations and low model stability). So procedure is:

$$
\begin{aligned}
relbias = (retro - base) / base \\
rho = mean(relbias)
\end{aligned}
$$

For SSB and F math approach will be

$$
\begin{aligned}
\rho_{SSB} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(B_{i} - \overline{B})}{\overline{B}} \\
\rho_{Fbar} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(F_{i} - \overline{F})}{\overline{F}} \\
\end{aligned}
$$

where i - year steps from terminal year (terminal-5, terminal-4, ... terminal).

```{r echo=FALSE}
retro.rho.ssbdata %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,1)) %>%
  knitr::kable(., caption = "Table 6.2.1. B retrospective values", row.names = TRUE)
```

```{r echo=FALSE}
retro.rho.fdata %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,3)) %>%
  knitr::kable(., caption = "Table 6.2.2. F retrospective values", row.names = TRUE)
```

**Final** Rho estimates: 

$$\rho_{SSB} = `r round(retro.rho.ssb,3)`$$
$$\rho_{Fbar} = `r round(retro.rho.f,3)`$$

## 7. Summary
There is summary plot from datalimited2 package displayed. This figures can be used directly to fishing regulation management.

```{r echo=FALSE}
plot_dlm(cmsy)
```

## Appendix list

### Appendix 1. Model reference points

```{r echo=FALSE}
print(cmsy$ref_pts)
```

### Appendix 2. Model estimates
```{r echo=FALSE}
print(cmsy$ref_ts)
```